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ABSTRACT 

Context. The image degradation produced by atmospheric turbulence and optical aberrations is usually alleviated using 
post-facto image reconstruction techniques, even when observing with adaptive optics systems. 

Aims. These techniques rely on the development of the wavefront using Zernike functions and the non-linear optimiza- 
tion of a certain metric. The resulting optimization procedure is computationally heavy. Our aim is to alleviate this 
computationally burden. 

Methods. To this aim, we generalize the recently developed extended Zernike-Nijboer theory to carry out the analytical 
integration of the Fresnel integral and present a natural basis set for the development of the point spread function in 
case the wavefront is described using Zernike functions. 

Results. We present a linear expansion of the point spread function in terms of analytic functions which, additionally, 
takes defocusing into account in a natural way. This expansion is used to develop a very fast phase-diversity recon- 
struction technique which is demonstrated through some applications. 

Conclusions. This suggest that the linear expansion of the point spread function can be applied to accelerate other 
reconstruction techniques in use presently and based on blind deconvolution. 

Key words. Techniques: image processing — Methods: analytical, numerical — Telescopes — Atmospheric effects 



1. Introduction 

Atmospheric turbulence degrades astronomical images by 
introducing aberrations in the wavefront. During the last 
years, functional adaptive optics systems have been devel- 
oped with the aim of partially cancelling these aberrations 
and generating a diffraction, limited image of any astronom- 
ical object (e.g., lBeckerslll993[ ). The development of such 
systems has been particularly challenging for solar tele- 
scopes because the reference object (the solar surface) has 
spatial structure. In spite of the great success of adaptive 
optics systems (Rimmele 2000; Scharmer et al. 2000), the 
limited number of degrees of freedom of the deformable mir- 
ror and the speed at which it has to work, especially in so- 
lar observations, impedes a full correction of the wavefront. 
Therefore, it has become customary in solar imaging to ap- 
ply post-facto reconstruction techniques to the observed im- 
ages to improve the wavefront correction. The advantage of 
these techniques is that the time limitation is much less re- 
strictive since they do not need to work in real-time. Hence, 
one can use lots of computation power to reconstruct the 
images. One of the methods considered first for the recon- 
struction of solar images w as based on speckle techniques 
(|von der Liihel 1 19931 Il994f ) which produced images of ex- 
traordinary quality. However, it suffers from two fundamen- 
tal problems. On the one hand, the number of images that 
one needs to acquire is large (although comparable with the 
number of images needed in complex blind deconvolution 
algorithms) because the missing information on the wave 
phases is statistically estimated from a succession of many 



short exposures. On the other hand, it relies on a theoret- 
ical description of the effect of the atmospheric turbulence 
and adaptive optics on the image ([Woger fc von der Liihel 
i2007t ), which might be inaccurate. However, recent efforts 
have shown that almost real-time processing (Denker et al.l 
12001 ^ and photometric precision (|Woger et al.l |2008[ ) can 
be achieved routinely. 

A complementary approach to speckle reconstruc- 
tion appeared with the systematic application of tech- 
niques based on p hase diversity (e.g. jLofdahl &: Scharmerl 



11994, 



Lofdahl et al. 1998|) , dev eloped some years before 
Gonsalves fc Chidlawl (|l979f ) and later extended by 
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iPaxman et al.l (|l992l ). These techniques rely on a maxi- 
mum likelihood simultaneous estimation of the wavefront 
and the original image (before being affected by atmo- 
spheric+telescope aberrations) using a reference image and 
a second one that contains exactly the same aberrations as 
the reference one plus a known artificially-induced aberra- 
tion. Because of the mechanical simplicity, a defocus typi- 
cally constitutes the added aberration. The main drawback 
of the phase diversity technique is the heavy computational 
effort needed to compute the maximum likelihood solution. 
Many Fourier transforms have to be used to estimate the 
point spread function (PSF) from the wavefront at the pupil 
plane. 

Even more computationally expensi ve is the blind 
deconvolution technique developed by Ivan Noort et al.l 
(j2005|). This technique, termed multi-object multi- frame 
blind deconvolution (MOMFBD), is able to reconstruct es- 
timates of the object and the wavefront from the observa- 
tion of different objects (same location in the Sun observed 
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at different wavelengths) with some constraints, which are 
essentially used to augment the amount of information in- 
troduced in the likelihood function. As a consequence, the 
estimation of the wavefront and the original image is less 
dominated by noise and hence more stable. In close analogy 
with phase diversity, this enhancement in stability comes 
at the price of a higher computational burden since many 
Fourier transforms have to be computed for the maximiza- 
tion of the likelihood. In conclusion, all these techniques 
suffer from time-consuming computations. As modern so- 
lar telescopes increase their flux of high quality data, ob- 
servers are condemned to struggle with the bottleneck of 
reconstruction algorith ms. 

During the 1940's, iNiiboerl (|l942f ) made it possible to 
estimate the PSF of an aberrated optical device in terms of 
the generalized complex pupil function. The resulting inte- 
gral expressions were very involved and remained impossi- 
ble to evaluate until recent years. What is presentl y known 
as the extended Nijbo er-Zernike theory (ENZ; iJanssenI 
I2QQ2I : iBraat et al.l l2QQ2"[ ) presents analytical expressions for 
the complex field at the focal volume that depend on the 
aberration through the coefficients of the Zernike expan- 
sion of the wavefront. The expressions are closed in the 
sense that they only depend on computable functions and 
the coefficients of the Zernike expansion. 

Our purpose in this paper is to extend the ENZ theory 
to arbitrary aberrations exploiting mathematical properties 
of the Zernike functions recently found. As a consequence, 
we find a general expression for the PSF in the focal volume 
as a linear combination of given basis functions. Starting 
from this expression, we propose to use it in a phase diver- 
sity reconstruction algorithm. Due to the analytical char- 
acter of the PSF, we are able to avoid the calculation of 
Fourier transforms during the iterative process, thus in- 
ducing an important gain in computational time. This shall 
be our fundamental result: that through this formalism we 
can drastically diminish the computational burden of re- 
construction algorithms and therefore the time needed to 
perform the reconstruction. We illustrate it using phase di- 
versity. But in principle any reconstruction algorithm can 
benefit from our formalism, since we only change the way 
in which the PSF is described. 

The outline of the paper is the following. Section [2] gives 
a description of the wavefront expansion in Zernike func- 
tions and how the field distribution can be analytically ob- 
tained for general aberrations. The expansion of the point 
spread function is presented in ^with an statistical anal- 
ysis of that expansion for atmospheric Kolmogorov turbu- 
lence in ^ Section [5] presents the application of the formal- 
ism to the fast restoration of images using phase-diversity. 
Finally, the conclusions are presented in ^ 

2. Wavefront expansion 

Let P(p, 0) represent the wavefront of a point source beam 
in a pupil plane, with p and the polar coordinates of 
the unit disk in that plane. The field distribution at an 
i mage p lane c an b e written in the Fresnel framework as 
(jB^rn L WolSll98Qh : 

U{r,ip;f) = - I I pdpd6>e^'^^'P(p,6>)e2^^^^^°^^^-^\ 
^ Jo Jo 

(1) 



which is, in general, a complex quantity. In that expres- 
sion we have defined a coordinate system (r, (p; f) over the 
image space, where (r, cp) are angular coordinates in the 
focal plane and / is the distance to that focal plane. We 
define the focal plane as the plane where the optical system 
would form the image of the object in the absence of any 
aberration. We implicitly assume an optical system with 
numerical aperture NA <C 1, without loss of generality for 
solar instruments. The (r, (p; f) coordinates can be put in 
relationship with the geometrical coordinates (X, F, Z) of 
the optical system by means of the following relationships 

2niNA) ^MNAl ^TriNAf 

r — J + tan(Z) = — , 

X 

(2) 

where A is the wavelength of the light. 
2.1. Wavefront expansion 

In order to solve the integral in Eq. ([T]), we have to pro- 
vide a useful description of the wavefront P(/), 0). In general 
the wavefront is also a complex quantity which is usually 
represented in polar form: the modulation in amplitude is 
represented by A{p^ 9) while the phase modulation is de- 
termined by <l>(p, 9). If the pupil plane is fully transmissive 
(or perfectly reflective) we can let ^) = 1 so that: 

P(p, 9) = A{p, 9)e'^^P^^^ = e^^^^'^) . (3) 

This is the case of a clear aperture telescope with the wave- 
front modified by atmospheric turbulence. For the sake of 
simplicity, we assume this to be the case, since it produces 
simpler mathematical expressions that can be related more 
transparently to previous results. An annular pupil would 
only require a change of the integration limits in Eq. ([T]), 
thus modifying accordingly the final result. 

It is customary (and for many reasons interesting) to 
express the phase modulatio n of t he wavefront in terms of 
the Zernike functions ( NqH [19761 ). The Zernike functions 
are labeled by two quantum numbers: a principal number 
n and an azimuthal number m fulfilling that \m\ < n and 
that n — \m\ is even. Their definition is: 

Z::{p,0) = R^{p)cos{m0) (4) 
Z-"^{p,0) = R^{p)sm{me) (5) 
= Rl{p\ (6) 

where R^{p) is a Zernike polynomial (e.g., iBorn fc Wolj 
19801 ) . In order to simplify mathematical manipulations, 
Nolll (|l976f ) introduced a useful single ordering index j , with 
perfect correspondence with pairs (n, m) so that higher val- 
ues of j represent higher aberrations with smaller probabil- 
ity amplitudes in a Kolmogorov turbulent atmosphere. We 
will use the two notations indistinctly. The reader should 
be aware that any pair (n,m) appearing expHcitly in an ex- 
pression has to be in agreement with the corresponding Noll 
index j. Since Zernike functions constitute a family of or- 
thogonal functions in the unit circle, any phase aberration 
can be written as a linear expansion: 

<^{p,e) = Y,a^z^{p,e) = Y.<^n{p,o). (7) 

j n,m 
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In a turbulent atmosphere of Kolmogorov type, the 
aberrations are characterized by coefficients aj that follow 
a multivariate Gaussian statistical distributions whose co- 
variance matrix can be obtained analytically (iNolllll976| ). It 
is interesting to point out that the covariance matrix of the 
coefficients is not diagonal, although it is possible to find 
modifications of the Zernike functions that diagonalize it. 
Returning to the expression of the exponential of the phase 
aberration, we expand it in power series: 



1 



(8) 



The radius of convergence of such expansion is infinite and 
that does not translate into an easy cutoff of the series. 
But using the decomposition in terms of Zernike functions 
of the phase aberration, we end up with: 



(9) 



For small aberrations {aj <C 1) the previous series can be 
safely cut at ffist order, s omething that has b een success- 
fully e xploited recently by ' JanssenI (|2QQ2f ) and iBraat et al.l 
(|2QQ2[ ) to empirically estimate point spread functions of 
microscope- type optical systems. Such approximation may 
be useful if one is interested in describing aberrations in- 
troduced by well adjusted optical systems, but it is hardly 
acceptable for aberrations introduced by atmospheric tur- 
bulence. Because of the complexity of the resulting expres- 
sion, there has been no interest in the past for this approach 
as a means to estimate the PSF of atmospheric seeing. 
The reason is that the series expansion of Eq. (|9]) gener- 
ates products ZjZk at second order and high-order prod- 
ucts Zj ■ ■ Zk at higher orders. In the absence of any clear 
rule to manipulate such products, the power expansion was 
rendered useless already at second or der. 

A recent mathematical result by iMatharl (|2QQ9f ]P1 has 
radically simplified the approach to that expansion and we 
use those results here for the ffist time to pursue the ana- 
lytical integration of the PS F of an?/ aberra tion introduced 
by atmospheric turbulence. iMatharl (|2QQ9f ) constructively 
demonstrated that the product of two Zernike functions 
can be written as a linear combination of Zernike functions 
(see Appendix [A|) . so that: 



Zj{p,e)Zj\p,e) = ^dkZk{p,9), 



(10) 



where the coefficients dk can be obtained by direct appli- 
cation of th e relations show n in Appendix [XI The product 
expansion of lMatharl (|2QQ9f ) can be applied recursively with- 
out difficulties, resulting in the following general result: 



(11) 



where the Ck coefficients can be calculated from the dk co- 
efficients of Eq. (p^Oj) . As a consequence, the expansion of 
Eq. can be written certainly as: 



1 



Y,PkZk{p^O), 



(12) 



where the /3k coefficients are complex in general. In other 
words, if the phase aberration can be expanded in Zernike 
functions, the wave front can also be expanded efficiently 
in these functions with complex coefficients /3k that can be 
inferred from the coefficients aj by a repetitive application 
of the coupling relation of Eq. ()A.2p and the coupling coef- 
ficient of Eq. ()A.3p . A direct consequence of Eq. ([12]) and 
Eq. ([9]) is that, for weak turbulence, the imaginary parts 
of the P coefficients are approximately equal to the Zernike 
coefficients of the expansion of the wavefront: 



Oik, 



(13) 



and strictly equal at ffist order. The weak turbulence condi- 
tion for that approximation often applies to our astronom- 
ical applications as demonstrated statistically in Section [H 
Through this approximation any method that allows us to 
estimate the value of the P coefficients from the PSF (like 
the phase diversity method shown below) directly gives us 
an estimation of the Zernike coefficients of the wavefront, 
at least to ffist order. 

2.2. The field distribution 

Having written the pupil function P(p, 6) as a linear com- 
bination of Zernike functions with complex coefficients 
the field distribution of Eq. ([1]) simplifies to a sum of inte- 
grals of the form: 

U{r,ip-J) = - / pdpdOe'^p'e^'^'P''''''^^-^^ 



JO 



+ y2/3j- I I " pdpdOe'^p' Zje^^'P'^^'^^^-^A) 

In the usual mathematical development of the ENZ the- 
ory we uncouple the two integrals above in radial and az- 
imuthal parts. The azimut hal part can be simplified us ing 
the Jacobi- Anger identity ([Abramowitz fc Stegun[[T972[ ): 

/'27r 

/ e^^iprcos{e-^)^ime^Q ^ 27rie^^^ (27rpr) , (15) 

where the Jmi'^Trpr) are Bessel functions of the ffist kind. 
Making use twice of the previous integral identity in Eq. 
([T4|) and reformulating the azimuthal part of the Zernike 
functions into a single e*^^ (that contains both the cos and 
sin functions into its real and complex parts for an unsigned 
m index), we end up with: 



/7(r, (/.;/) = 2 f pdpe^^^" M2^pr) 
Jo 

+ 2^z-e^-^ft-F-(r,/). 



(16) 



The function V^{r^f) represents the following (complex) 
radial integral: 

/) = /' pdpe'fp"R\r^{p)J^{27Tpr), (17) 

JO 



k=2 



which fulfills that: 



See also http://arxiv.org/abs/0809.2368 



[c^(r,/)]^ = Kr(^,-/)- 



(18) 
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Fig. 1. Example of unique basis functions for representing the PSF including terms up to the Noll index j = 7. These 
functions are calculated for a telescope of 100 cm diameter at a wavelength of 5250 A, a typical situation in present day 
telescopes. The spatial dimensions are in units of pixels, that we choose of size 0.055", also typical in present observing 
conditions. The upper panel shows the results for the focused image while the lower panel shows the results for the 
defocused image, where the defocus is chosen as shown in ^ 



If we are interested in the field distribution at the plane has an easy solution: 
where the image would form in the absence of any optical 

aberration (/ = 0), the radial integral becomes real and V'^(r) = ( 1) "^"j"^' ^n+i(27rr) 
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Prefiguring its apphcation in phase diversity, the integral of 
Eq. (pT|) can be numerically calculated for arbitrary values 
of / using the ra pidly convergent series expansion developed 
bv lBraat et all (2QQ2h . 

The final expression for the field distribution is given 

by: 



Uir,ip;f) = 2V°ir,f) + 2j2 



I e 



^PjV:^ir,f). (20) 



Notice that, when no aberrations are present, the in- focus 
field distribution (at / = 0) is proportional to: 



Ji(27rr) 
27rr 



(21) 



which is the function giving rise to the well known Airy 
spot for a circular non- aberrated aperture. At / 7^ 0, the 
field distribution is proportion al to Vn(r,f ) whic h, using 
the rapidly convergent series of iBraat et al I (|2QQ2f ). can be 
written exactly as: 



00 

^=1 



(-2i/)' 



-1 J;(27rr) 
(27rr)' 



(22) 



Summarizing, the product expansion of'Mathar' (2009) 
has allowed us to write the field distribution in the image 
space as a linear combination of functions V^{r^ f) defined 
above, and of which the (in-focus or out-of-focus) Airy spot 
is just the first contribution. This can be c onsidered to be 
a ge neralization of the results presented by IJansse ^ (120021) 
and lBraat et al.l (|2002[ ) for the case of small aberrations. 



The functions appearing in the linear expansion can be cal- 
culated as: 



K,;:^ir,^) = (-1)™ C-'™ (^)Kr(r)y„T (r) 

i^C;™^^^^') = -(-i)™'5--™'(^)TC(0v;'^'(r)- 



(26) 



Note that the azimuthal and radial contribution are sepa- 
rated. The angular dependence is modulated by the follow- 
ing functions: 



COS 



TT 



(m + m') + (m — m')if 



— (m + m') + (m — m')Lp 



(27) 



The coefficients Cjk and Djk in the linear expansion 
of the PSF are related to the real and imaginary parts of 
products of P coefficients: 



Cjk 
Djk 



-Dkj = ejklm{f3jf3l) 



(28) 



where Ckj = 1 — Skj/2. As an example. Fig. [T] shows the 
unique basis functions up to Noll index j = 7. The rest 
of functions with different combinations of indices can be 
related to these by trivial relations like those shown in 
Appendix [Bl Note that the terms j = k of the summa- 
tion contain only the contribution of the Cjk coefficients 
and are azimuth-independent, since: 



n,n 



0. 



(29) 



3. The point spread function 

The point spread function or intensity distribution is ob- 
tained from the complex field distribution as: 



PSF(r,^; f) = Uir,ip; f)U^{r,^; /). 



(23) 



For the sake of simplicity in the presentation, we separate 
the analysis of the in-focus and the out-of-focus PSFs. 

3.1. In-focus PSF 

Using Eq. (|2Q|) with the first term included in the sum- 
mation as j = 0, the PSF results in the following double 
summation: 



PSF(r,^) = 4V 



'^'^Pj{-ir'e-'^'^(3lv;^{r)V;:^'{r), 

(24) 

where the (n, m) pair is associated with the Noll index j 
while the {n\ m') pair is attached to the Noll index k. The 
previous double summation can be separated into two sum- 
mations: one containing terms for which j = k and one 
containing terms j ^ k. After tedious but straightforward 
algebra, the final expression for the point spread function 
is: 

PSF(r, ^) = 8 ^ [Q.^n^;T^>, ^) + ^.^Cn^V, ^)] .(25) 

k<j 



Because n — \m\ is even in the definition of the Zernike 
functions, it can be verified that only the terms with j = k 
in Eq. ([25]) contribute to the total area of the PSF, so that: 



/ / rdrd(pPSF(r, (/:))= V 

^0 Jo • n 



1/3. 



7r(n + 1) 



(30) 



where the value of n is associated to the Noll index of sum- 
mation. As a result, the Strehl ratio defined as the ratio 
between the peak intensity at the origin of the PSF and 
the equivalent Airy function simplifies to: 



(n + 1) 



(31) 



3.2. Out-of-focus PSF 

Following the same approach as in the previous section, we 
can also write an expression for the PSF at any arbitrary 
focal position around the focal point. Since in this case 
the V^{r^f) functions are complex in general, the basis 
functions are slightly more elaborate: 

PSF(r, ^; /) = 8 ^ [Q^^^h'^V, /) + D,kY:^f (r, ^; /) 



k<j 



(32) 
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Fig. 2. Histograms of the distribution of the wavefront expansion coefficients a for an atmosphere with Kolmogorov-type 
turbulence (upper row). The real and imaginary part of the corresponding expansion coefficients /3 of the field distribution 
are shown in the middle and lower rows, respectively. The simulations are carried out for a telescope of I) = 90 cm and 
a Fried parameter tq = 7 cm. 



where 

^r;™'(^^;/) = (-ir 



Cn"(^^;/) = (-ir 

(33) 

The auxiliary functions appearing in the previous formulae 
are given by: 

/) = Im(Kr)Im(yj') + Re(yr)Re(yy ) 

^Z'^^'ir; f) = Re(TC)Im(Fy ) - Im(TC)Re(TC' ). (34) 

Obviously, when the previous expressions are evaluated 
at the focus (/ = 0), because the functions are real, 

'^rTn^ = and ^^^7' = y^y^ , thus recovering the ex- 
pressions of Eq. ([26]) . Likewise, for the terms of the sum- 
mation in Eq. (|32]) for which j = /c, we have: 



^; /) = (n /) = (r; /)[TC(n /)]^ 

>rr(^'<^;/) = o- (35) 

4. Statistical properties of the /3 coefficients 

In a Kolmogorov turbulent atmosphere, the coefficients otj 
of the Zernike expansion of the wavefront follow a multivari- 
ate Gaussian distribution with a non-diagonal covariance 



mat rix determined by iNollI (|1976[ ) , IWang fc MarkevI (|l978f ) 
and Roddier (1990). The amplitude of the elements of the 
covariance matrix is determined by the Fried parameter tq. 
As a function of this parameter, we simulate random sam- 
ples Qffc from such multivariate Gaussian and use the expres- 
sions of Section [2] and the formulae developed by Ma tharl 
(2009) presented in Appendix El to calculate the equiva- 
lent coefficients. The corresponding statistical analysis 
is presented in Fig. [2] for aberrations of order higher than 
two, excluding tip- tilt. The upper row shows that the a/^ 
coefficients are Gaussian distributed, with a variance de- 
pending on the order of the aberration. The middle and 
lower rows show the histograms for the real and imaginary 
part of the coefficients, respectively. 

The distribution of the imaginary part is close to 
Gaussian and very similar to that of the wavefront ex- 
pansion coefficients. As already mentioned, this is a con- 
sequence of the fact that, to first order, the imaginary part 
of the coefficients equals the a/c coefficients. Higher- 
order corrections to the imaginary parts occur only at odd 
orders, so that the first correction takes place already at 
third-order, which often induces a relatively small change 
for turbulence not too strong. The Pearson linear correla- 
tion coefficient between a/c and Im(/3/c) is larger than 0.95 
for practically all orders, indicating that Im(/3/c) represents 
a good estimation of the wavefront expansion coefficients. 
The real parts present distributions with a smaller variance 
and, in some cases, a large kurtosis. This is a consequence of 
the fact that the real part of the coefficients is modified 
only at second-order often representing a small correction. 
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5. Phase diversity 

A fundamental advantage of the analytical expression for 
the PSF just presented is that it allows to introduce a 
defocus in the PSF while maintaining unchanged the Pj 
coefficients that characterize the aberrations. As a con- 
sequence, it is natural to propose such a functional form 
for doing ima^e reconstruction first through phase diver- 
sity techniques (Gonsalves fc C hidlawlllQTQclPaxman et all 
li£92; Lofdahl & Scharmer 1994b^ In order to fix naming: 
conventions, we first summarize the basics of the standard 
approach to phase diversity. 



5.1. Standard approach 

The idea behind phase diversity is to simultaneously ac- 
quire a focused image (I^o) and another one with exactly 
the same aberrations plus a known defocus (Dk). Under 
the standard image formation paradigm, the observed im- 
age after degradation by the atmosphere and the optical 
devices of the telescope can be written as: 



Do = F^Po^No 



(36) 



where Pq and Pk are the PSFs of the atmosphere+telescope 
in the focused and defocused images, respectively. The 
quantities A^o and Nk are noise contributions in the image 
formation, each one characterized with a variance ctq and 
(7^, respectively. Since both images are obtained simultane- 
ously, the underlying object F is the same in both cases. 

In the standa rd phase div ersity method (see , 

e.g.. iGonsalves fc Chidlaw 119791 : iPaxman et all 119921 : 
iLof dahl fc Scharmer 1994b), the phase aberration in the 
focused image is expanded in a Zernike basis according to 
Eq. ([710. The PSFs are obtained from this phase aberration 
by calculating: 

Po = \FT-' {A{p,0)exp[i^p,0)]}\' 

Pk = |FT-^{A(p,^)exp[^$(p,^)+^(5,(p,^)]}|^ (37) 

where the defocus is introduced by adding a term of the 
form 5d = afZ/^^p^O) and FT""*^ stands for the inverse 
Fourier transform. The election of the defocusing distance 
is essential for obtaining an efficient and successful phase 
diversity method. It is usually chosen to represent one wave 
retardance at the wavelength of the observation. 

The coefficients of the expansion in Zernike functions of 
the phase aberration are obtained by calculating the coef- 
ficients aj of the expa nsion of Eq. tha t minimize the 
following error metric ([Paxman et al jTl992f ) : 



|Po|2+7|A|2)l/2' 



(38) 



where D^^ Dk^ Po and Pk are the Fourier transforms of I^o, 
Dk, Pq and P/c, respectively. The summation is carried out 
over the full Fourier domain defined by the frequencies u 
and while 7 = (JQ/al. Note that the dependence of the 
metric on the aj coefficients is highly non-linear due to the 



^ Some advantage is gained by using instead the Karhunen- 
Loeve expansion. The reason is that these functions are empiri- 
cally built to explain the largest amount of turbulence variance 
with the fewer number of functions. 



exponential function, the Fourier transforms and the mod- 
ulus operation that appears in the definition of the PSF. 
Note also that each evaluation of the metric and its gradient 
with respect to the aj coefficients needs to carry out sev- 
eral Fourier transforms, which are of order 0{N log N). At 
a final step, once the PSF (and the ensuing Fourier trans- 
form) are known, the reconstructed image is estimated as 
the inverse Fourier transform of (Pax man et al 1992) : 



DoP^^jDkPl 
|PoP+7|AP ' 



(39) 



The maximum- lik e lihood solution of the metric of 
IGonsalves fc Chidlawl (|l979f ) is not affected by the pres- 
ence of additive Gaussian noise but it strongly affects the 
deconvolution process given bv Eq. (|39|) since noise is un- 
Hmitedly amplified. We follow iLofdahlfc Scharmerl (|l994b[) 
and apply a filter H in the Fourier domain to the Fourier 
transform of the observed images with the aim of filtering 
a la rge part of the n o ise. Th e same procedures applied by 
Lofd ahl fc Scharmed (|l994b[ ) for filtering remammg noise 
peaks at high frequency are also applied. 

5.2. The new approach 

The crucial advantage of the phase diversity algorithm we 
propose is that the PSF of the focused and the defocused 
images can be written as linear combination of known func- 
tions. Because of that, the previous error metric can be 
minimized without performing a single Fourier transform 
during the iterative scheme. After Fourier- transforming the 
observed images and the basis functions X and F, it is a 
matter of finding the values of the /3j complex quantities 
that minimize the metric L. Another advantage is that the 
non-linear dependence of the metric function on the Pj co- 
efficients is less pronounced, thus helping the optimization 
routine to find the minimum. 

In order to find the values of the coefficients Pj that 
minimize t he metric L, we use a scaled conjugate gradient 
algorithm (|Andrei|[2QQ 8) that performs extremely well and 
fast for complex problems. To this end, one needs the gradi- 
ent of L with respect to the Pj coefficients. After a tedious 
but straightforward algebra, the gradient can be found to 
be: 



dL 



dp^ 



- = 2yRe 

r,i Z-^ 



\ "i Pi J 



2o„ I £,t , pt 



\ Pi Pi 



(40) 



where Q = (|PoP + l\Pk?)-^'^ and: 

k=0 

E (/5fetTV' 0) + /5i^™;'^V' 0)) 

k=0 
k=0 

E {PlK^'i^^ /) + l^iKn^'ir, <P; /)) . (41) 



dPo 

P\ 
dpQ 

P\ 

dPk 
PI 
dPk 



k=0 
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The superindices r and i indicate the real and imaginary 
part of the coefficient /3i , respectively. One should be aware 
that the values of n and m in the previous equation are 
those associated to the Noll index while n' and m' are 
those associated to the index k of the summation. Finally, 
the functions X and Y are the Fourier transforms of the 
functions defined in Eq. (|33l) . 

The computer code is written in Fortran 90 and runs in 
parallel environments using the Message Passing Interface 
library. Prior to the phase diversity restoration procedure, 
there is an initial phase in which the basis functions of Eqs. 
(|26l) and (|33]) are precalculated and their Fourier transforms 
are saved for later use. Since this has to be run only once for 
each telescope and size of the isoplanatic patch, it can be 
done with the required precision at virtually no extra com- 
puting time. In a second phase, in which the actual phase 
diversity restoration process is carried out, these Fourier 
transforms of the basis functions are quickly read and used 
in all subsequent calculations. This greatly accelerates the 
computing time and only two Fourier transforms are needed 
per isoplanatic patch. The parallel implementation scales 
very well with the number of processors. The field-of-view 
is divided into isoplanatic patches of a given size and the 
reconstruction of each patch is carried out by one of the 
available processors. In principle, if A^proc are at hand, the 
computing time is reduced by a factor l/A/'proc because the 
contribution of all input/output tasks is negligible with re- 
spect to the actual iterative process to minimize the phase 
diversity merit function. 

5.3. Illustrative examples 

We present two representative examples of the success of 
the implementation of our analytical representation of the 
PSF for carrying out phase-diversity image reconstruction. 
We stress that our purpose is not to present amazing re- 
constructions but to demonstrate that reconstruction can 
be made to the same level of quality of the standard pro- 
cedures only much faster. 

5.3.1. Blue images 

In our first example the quiet Sun internetwork was 
observed with the S wedish 1-m Solar Telescope (SST; 
IScharmer et"all |2QQ2[ ) during 2 007 September 29 during a 
period of very good seeing (see Bonet et al. 2008, for a de- 
scription of the observations). The pixel size is 0.034^' and 
the observations are carried out in the continuum close to 
the Ca H line at 3953 A, providing a theoretical diffraction 
limit of 0.1 '^ All these parameters were chosen to enhance 
granulation contrast. With an isoplanatic patch of 128 x 128 
pixels {^4.5" X 4.5'^) , an image of 1000 x 1000 pixels contains 
approximately 60 isoplanatic patches. Our phase diversity 
restoration algorithm takes around 5 s per patch when in- 
cluding 20 terms in the summation of Eq. (p!2]) that de- 
scribes the wave front. Therefore, the total computing time 
is ~5 min. When run in a standard present-day desktop 
computer with 4 processors, every image can be fully re- 
stored in a timescale of 1 minute. Fig. [3] presents the results 
of the reconstruction (right panel) and the original focused 
image (left panel). The region zoomed illustrates the spatial 
frequencies enhanced by the phase- diversity algorithm. 



5.3.2. Red images 

In this second example, we analyze observations of an active 
region carried out during 2002 May 15 with the THEMIS 
telescope (Rayrole & Mein 1993) at a wavelength of 850 
nm. The pixel size is 0.075'^ and, due to the large wave- 
length, the diffraction limit is 0.24^', limiting the contrast 
of the granulation, both because of the wavelength and the 
intrinsically reduced contrast in the infrared. This dataset 
has been also use by Criscuol i et al.l (2005) to present the 
first results of phase-diversity applied to THEMIS data. 
The images were restored with patches of 100x100 pixels, 
equivalent to 7.5" x7.5" . Patches were finally mosaicked to 
build the final image shown in Fig. SI were only the part 
inside the rectangle has been reconstructed, using 20 terms 
in the summation of Eq. (p!2]) that describes the wavefront. 
Visual inspection clearly indicates an improvement in the 
image quality. If desired, better visual and quantitative re- 
sults can be obtained by adding multi-image inforni ation 
to the reconstruction algorithm (Criscuoli et al.ll2QQ5l ). 

Computing time is what makes the difference, since each 
patch, as in the previous case, only took 5 sec in a standard 
desktop computer. With such requirements, image recon- 
struction can be easily implemented as on line task in the 
data pipelines of the telescopes themselves. 

6. Conclusions 

The use of the recently developed Extended Nijboer- 
Zernike theory together with some mathematical results on 
the multiplication of Zernike polynomials has allowed us 
to rewrite the image formation integrals with an aberrated 
wavefront in terms of linear combinations of analytic func- 
tions. Generalizing the usual ENZ theory, we are able to do 
so a priori independently of the amplitude of the aberra- 
tions. With such mathematical tool we are able to rewrite 
the techniques of post-facto image reconstruction taking 
advantage primarily of the fact that the Fourier transforms 
of those analytic functions can be precomputed once and for 
all. Building different PSFs in the minimization process at 
the core of these reconstruction algorithms does not require 
to recompute any other Fourier transform but just modify 
the scalar, although complex, coefficients multiplying the 
functions. 

The gain in speed is enormous. We illustrate it through 
two observations readied for phase diversity. Image recon- 
struction with quality identical to the standard algorithms 
is performed in question of seconds per patch, and one to 
five minutes for wide-field images, in standard desktop com- 
puters. If image reconstruction is a must in present and 
future solar observations, the bottleneck of computation- 
ally expensive and time consuming algorithms was almost 
a showstopper for future instruments. The present improve- 
ment through the use of analytical PSFs solves the core of 
that problem and suggests that image reconstruction can 
even be implemented as a routine on-line procedure on the 
data pipelines of the telescope instruments themselves. 

Since we have only modified the description of the PSF 
with respect to previous approaches, all common recon- 
struction schemes (with any desired complexity), and not 
just phase diversity, can be extended to use our formalism. 
Among t hem, we find me thods like multi-image phase diver- 
sity (e.g.. IPaxman et al.lll 992) or multi-obj ect multi- frame 
blind deconvolution (|van Noort et al.ll200'5[ ) that use many 
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Fig. 3. Reconstruction example at 395 nm. The lower image shows the in- focus image while the right image presents the 
reconstructed image using only one pair of in- focus and defocused images. The contrast of the original image is 9.6% 
while the reconstructed one presents a contrast of 11.9%. Data has been acquired with the SST telescope. 




Fig. 4. Example of the phase-diversity reconstruction at 
850 nm. We represent the focused image (left) and the re- 
constructed sub- image indicated with a box (right). Images 
have been acquired with THEMIS. 



images to estimate information about frequencies that has 
been destroyed by the presence of noise. In essence, all these 
schemes require writing an error metric (equivalently, a like- 
lihood function) like Eq. (|38|) that takes into account the 
presence of the additional information. The larger amount 
of information helps regularize the problem and reduces the 
influence of noise. All of them are penalized by the comput- 
ing time involved and all of them can potentially take proflt 
of the analytical approach we present to reduce the compu- 
tational burden. The flnal goal is not to slow down the flows 



of images from instruments because of the lack of capability 
to handle the image reconstruction problem. 

Beyond those known techniques, we also want to point 
out that it is possible to introduce regularization by us- 
ing a fully Bayesian approach in which prior information 
is introduced in the problem, eliminating a priori the noise 
nuisance in the minimization of the metric. 
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Appendix A: Product of Zernike functions 

The breakthrough brought over bv lMatharl (|2QQ9 [) has been 
to point that, in the product of two Zernike functions, only 
the following three situations could arise: 



2 cos(mi^) cos(m2^) 
2 sin(mi^) cos(m2^) 
2 sin(mi6>) sm{m20) 



cos[(mi 
sin[(mi 



m2)0] 
m2)e] - 
cos[(mi — 1712)0] 



- cos[(mi + m2)0] 

- sin[(mi + m2)0] 

- cos [(mi -^ii^B\) 



In other words, the result is always the sum of trigonometric 
functions depending on two azimuthal numbers mi ± m2. 
This is indeed the seed of two new Zernike functions with 
ms = mi ± m2 if we succeed in writing the product of the 
radial polynomials with identical ms azimuthal value. That 
is, we seek a product expansion for the radial polynomials 
of the form: 



ni+n2 

E 

n3=m3 



i/ni, mi, 77,2, m2,n3,m3-^ 1^77,3 • 



(A.2) 



Every term of the series above has the same azimuthal 
number ms, coincident with the ms of the trigonomet- 
ric product and can therefo re be comb ined with it to 
build a Zernike function Zg^^ iMatharl (|2QQ9. ) demonstrates 
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that this is indeed feasible and provides the value of the 
5'ni,mi,n2,m2,n3,m3 Coefficient in the product expansion as 

^ni,mi,n2,m2,n3,m3 ~ 2(713 H~ 1) 
— ai —a2 —as 

X EEE 



1 



Si=0 S2=0 S3=0 

3 



^1 + ^2 + ^3 + 2(1 - Si - S2 - Ss) 



X n(-i) 



(A.3) 



where a = — (n — m)/2. These coefficients play the same 
role of the Clebsh-Gordan recoupling coefficients for the 
product of two spherical harmonics. 



Appendix B: Some properties of Xj;^'? and Y"^^ 
and analytical Strehl ratio 

Under the interchange of the {n^m) and in' ^m') pair of 
indices, the following expressions hold: 

c;r(^;/) = 0'^'(^;/) 

C^''^((^) = C^'^'(-(^) 

5^''^((^) = 5^'^'(-V^)- (B.l) 

As a consequence, the functions X^^'J^ and fulfill 
the following properties: 
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>;^'r(r,<P;/) = -(-l)™'-"X„":„T'(r,-^;/). (B.2) 



Concerning the Strehl ratio it is possible to demonstrate 
that only the Airy part of the PSF (either in-focus or out- 
of-focus) contributes to the central peak of the PSF. The 
reason is that, evaluating Eq. (p!?]) at r = gives: 



so that 



SnoSmo^ [sin / + i(l - cos/)] 



(B.3) 



0, 



1 — COS / 

2/2 



(B.4) 



which gives, for the / = 0, a contribution of 1/4 at the 
center of the PSF. 
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